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In numerical simulations, spontaneously broken symmetry is often detected by computing two- 
point correlation functions of the appropriate local order parameter. This approach, however, com- 
putes the square of the local order parameter, and so when it is small, very large system sizes at high 
precisions are required to obtain reliable results. Alternatively, one can pin the order by introducing 
a local symmetry breaking field, and then measure the induced local order parameter infinitely far 
from the pinning center. The method is tested here at length for the Hubbard model on honey- 
comb lattice, within the realm of the projective auxiliary field quantum Monte Carlo algorithm. 
With our enhanced resolution we find a direct and continuous quantum phase transition between 
the semi-metallic and the insulating antiferromagnetic states with increase of the interaction. The 
single particle gap in units of the Hubbard U tracks the staggered magnetization. An excellent data 
collapse is obtained by finite size scaling, with the values of the critical exponents in accord with 
the Gross-Neveu universality class of the transition. 

PACS numbers: 71.27.+a,71.10.Fd,71.30.+h,73.21.Ac,75.70.Cn 



I. INTRODUCTION 

Detecting spontaneous symmetry broken phases in nu- 
merical simulations often relies on the measure of corre- 
lation function. For instance, the magnetically ordered 
phase is characterized by long ranged spin-spin correla- 
tions, whereas the superconducting state exhibits long 
ranged pair correlations in the appropriate symmetry 
channel. A fundamental caveat with such an approach 
is that one measures the square of the order parameter. 
If the later is small, the quantity one attempts to obtain 
by extrapolating numerical data to the thermodynamic 
limit is quadratically smaller. As a consequence, very 
large system sizes at high precision are required in addi- 
tion to an appropriate finite size extrapolation formula. 

The aim of this article is two-fold. We will document 
a simple and very efficient alternative method to detect 
magnetically ordered phases in SU(2) invariant Hubbard 
type models in the realm of projective quantum Monte 
Carlo methods. With the enhanced resolution, we will re- 
visit the semi-metal to insulator transition on graphene's 
honeycomb lattice, which has recently been under con- 
siderable debate [TJI2]- 

Form the technical point of view, our approach is very 
similar in spirit to an approach considered in Ref. El 
By introducing a local magnetic field at say the origin, 
we explicitly break the SU(2) spin symmetry. In the 
presence of long range order and in the thermodynamic 
limit, any field will pin the order along the direction of 
the external field. Thereby, order can be detected by 
computing directly the magnetization infinitely far from 
the pinning field. The upside of such an approach is that 
one measures directly the order parameter rather than 
its square. This amounts to evaluating a single particle 



quantity which is often much more stable than correlation 
functions. The downsides are three-fold. One explicitly 
breaks SU(2) spin symmetry such that spin sectors mix 
and it becomes computationally more expensive to reach 
the ground state. Since the computational cost scales 
linearly with the projection parameter, this problem is 
tractable. The second difficulty lies in the ordering of 
limits. To obtain results which are independent on the 
magnitude of the pinning field, it is important to first 
take the thermodynamic limit and then the limit of infi- 
nite distance from the pinning field. In a practical imple- 
mentation, this ordering of limits has as a consequence 
some leftover dependence of the magnetization on the 
magnitude of the pinning field. This is particularly visi- 
ble when the pinning field is small. The final drawback 
is that it is not always possible to introduce a pinning 
field without generating a negative sign problem. For 
instance, in the Kane-Mele Hubbard model |U [5], the 
spin order lies in the x-y plane. Adding a magnetic field 
along this quantization axis introduces a sign problem. 
On the other hand, the method is applicable to SU(N) 
symmetric Hubbard-Heisenberg models [6 . 

The organization and main results of the article are the 
following. We will focus on the Hubbard model on honey- 
comb lattice at the filling one half, for which the presence 
of an intermediate spin-liquid phase has been controver- 
sial [TJ |2]. After introducing and testing the approach 
in the next section, we will provide a phase diagram of 
the Hubbard model in Sec. |HI| The data points to the 
fact that the staggered moment follows rather precisely 
the single particle gap, when the latter is measured in 
the natural units of the Hubbard U, suggesting a direct 
quantum phase transition between the semi-metallic and 
the insulating antiferromagnetic phases. Furthermore, an 



excellent finite size scaling of the data for both the stag- 
gered magnetization and the single particle gap is found 
by assuming the values of the critical exponents (3 = 0.79 
and v = 0.88. These are the values found in the first or- 
der expansion for the Gross-Neveu- Yukawa field theory of 
this quantum phase transition [7], around its upper crit- 
ical (spatial) dimension of three. Altogether, the data 
strongly supports the existence of a single quantum criti- 
cal point separating the semi-metallic and the insulating 
antiferromagnetic phases of the Hubbard model, with the 
quantum criticality belonging to the Gross-Neveu univer- 
sality class [Bi- 



ll. MODEL AND METHOD 

As in Ref. pQ, we will consider the half-filled Hubbard 
model on the Honeycomb lattice 

Hm = -* £ c{ a c ha + C/^Kt - V2) («u - V 2 ) • 

(1) 

The hopping is restricted to nearest neighbors such that 
the bipartite nature of the lattice allows us to avoid the 
negative sign problem. 

Generically, to detect anti-ferromagnetic ordering we 
compute spin-spin correlations: 
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Here N — 2L 2 corresponds to the number of orbitals, 
and L is the linear length of the lattice. A finite value of 
to signalizes long range order and is equivalent to spon- 
taneous symmetry breaking. In particular, including a 
magnetic field term with appropriate Fourier component, 



H h = hJ2e lQi S*, 
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FIG. 1: Comparison of the pinning field and correlation 
function approaches to determine the staggered moment. The 
data sets at ho = corresponds to the correlations functions, 
and values of 9t = 40 are sufficient to converge to the ground 
state. For non-vanishing pinning fields projection parameters 
6t = 320 are required to guarantee convergence. We have used 
Art = 0.1 which for the symmetric Trotter decomposition 
yields converged results within our numerical accuracy. At 
U/t — 4 (b) comparison with results of Ref. [2] shows excellent 
agreement. Lines corresponds to least square fits to the form 
a + b/L + c/L 2 . 
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The ordering of limits is crucial. One first has to take 
the thermodynamic limit to allow for the collapse of An- 
derson's tower of states, and then the limit of vanishing 
magnetic field h. 

It is more convenient to consider a local field since, as 
we will see below, this lifts the burden of taking the limit 
h — > numerically. The local pinning field is given by 
the term 



Hi„r — hnS, 



o^o 



(5) 



in the Hamiltonian. Using the representation 5i t o — 
jj2 X) q e * q °f t ne Kronecker symbol shows that each 



that taking the thermodynamic limit is equivalent to tak- 
ing the amplitude of the relevant Fourier component to 
zero. With the local field construction, the appropriate 
ordering of limits for anLxi lattice reads: 
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That is, one first has to take the thermodynamic limit 
- again to guarantee the collapse of the tower of states in 
the presence of long range order - and only then can one 
take the distance from the pinning center to infinity [15] . 
In other words, the distance from the pinning center sets 
an energy scale which has to be larger than the finite size 
spin gap. As an efficient estimator for the evaluation of 



the ordered moment we thus propose: 

m = r lim j2J2 elQ ' i ^ H ^+ H ^ 



(7) 



We have tested the above approach for the Hubbard 
model on the honeycomb lattice. Ground state calcula- 
tions were carried out with the projective auxiliary field 
Quantum Monte Carlo (QMC) algorithm which is based 
on the equation: 
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Here 9 is a projection parameter, and the trial wave func- 
tion is required to be non-orthogonal to the ground state. 
For H — H t (j + Hi oc the inclusion of the magnetic field 
does not generate a negative sign problem. We have cho- 
sen the trial wave function to be the ground state of 
the non-interacting Hamiltonian Ht = H t + Hi oc in the 
S z = sector. The implementation of the algorithm fol- 
lows closely Refs. [TJ |S| The major difference is the use 
of a symmetric Trotter breakup which ensures hermitic- 
ity of the imaginary time propagator for any value of the 
time discretization At. It also leads to smaller system- 
atic errors. 

Figure fTJ (a) plots the local moment at U/t = 5 us- 
ing different methods. In this case, magnetic ordering is 
robust such that various approaches can be compared. 
The data set at ho = corresponds to the correlation 
functions of Eq. [2j For this set of runs we used a spin- 
singlct trial wave function and the projection parameter 
8t = 40 suffices to guarantee convergence to the ground 
state. This quick convergence stems from the fact that 
the trial wave function is orthogonal to the low lying spin 
excitations [10 . The runs at finite values of the pinning 
field correspond to the quantity of Eq. W\ In the pres- 
ence of a finite pinning field, SU(2) spin symmetry is 
broken and the trial wave function overlaps with all spin 
sectors. Consequently, a large value of the projection pa- 
rameter Ot = 320 is required to guarantee convergence to 
the ground state within the quoted accuracy. Note that 
the CPU time scales only linearly with the projection pa- 
rameter, so that such large projection parameters are still 
numerically tractable. It is also worth pointing out that 
the observable of Eq. [7] corresponds to a single particle 
quantity and shows very little fluctuations. As evident 
in Fig. lija), finite size effects are strongly dependent on 
the specific choice of the pinning field. Nevertheless, con- 
vergence to values consistent with the generic approach 
based on Eq. [2] is obtained for relatively large values of 
the pinning field. If the pinning field is chosen too small, 
larger lattices are apparently required to insure that the 
finite size spin gap is smaller than the energy scale set 
by the pinning field. This expectation is confirmed by 
the data, which shows a systematic upturn as a function 
of the system size for smaller values of the pinning field. 
Such a non-monotonous finite size behavior complicates 
a finite size scaling analysis. For this reason, we propose 
to use a relatively large value of the pinning field. 



At U/t = 4, the local moment is smaller, and hard to 
detect. The ho = data set of Fig. [l](b) compares our 
results for the correlation function of Eq. [2] to those of 
Ref. [5J As apparent, the agreement up to our largest 
lattice size, L — 18, is remarkable. Without the largest 
lattice sizes, L = 24 and L = 36, extrapolation to the 
thermodynamic limit is hard due to the downward turn 
present in the finite size results. The data sets stemming 
from the pinning field approach provide an alternative 
perspective, and, on the whole, confirm the result of Ref. 
2, For the considered field range there is considerable 
scatter in the finite size results, but nevertheless the ex- 
trapolation to the thermodynamic limit seems to be field 
independent, as expected from the above considerations. 

We conclude this section by mentioning that we have 
tested the approach for the non-interacting case and the 
method successfully demonstrates the absence of long 
range magnetic order. Hence both for critical states as 
well as for magnetically ordered phases the pinning field 
approach does provide an efficient tool. Further testing of 
this approach for Hcisenbcrg bilayers is presently under 
progress [TB] . 



III. PHASE DIAGRAM OF THE HUBBARD 
MODEL ON HONEYCOMB LATTICE. 

We have used the above approach to revisit the mag- 
netic phase diagram of the Hubbard model on the hon- 
eycomb lattice. At weak couplings the model is known 
to have a stable semi-metallic state. In the strong cou- 
pling limit, and due to the absence of frustration, an 
anti- ferromagnetic Mott insulator is present. The nature 
of the transition between these two states has been stud- 
ied in the past, [HI Q3J [H] and is presently controversial 
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A. Single particle gap. 

To pin down the coupling strength beyond which the 
single particle gap opens, we have repeated calcula- 
tions for the time displaced single particle imaginary 
time Green function at the nodal point: G(K,t) — 

E ct (4,» c k,> = 0)). As evident in Fig. B and 
with the symmetric Trotter decomposition, the Trotter 
systematic error is negligible within our accuracy. Fit- 
ting the data to an exponential form allows us to extract 
the single particle gap, which we plot as a function of 
system size in Fig. [21 Assuming a polynomial form for 
the extrapolation to the thermodynamic limit, we find a 
small but finite single particle gap for U/t > 3.7. This 
finding is particularly interesting when compared to the 
results of Sorella et al [5] which at U/t = 3.8 point to 
the absence of long range magnetic order. This could be 
taken as an indication for a possible intermediate phase. 
However, the analysis that we present below suggests a 
different interpretation. 






<r 




FIG. 2: Single particle gap. (a) The raw data at U/t = 
3.8. As apparent, the systematic error stemming from the 
finite Trotter step is negligible within our accuracy. The data 
supports a large imaginary time range consistent with a single 
exponential decay. Lines are least square fits of the tail of the 
imaginary time Green function to the form Ze~ spT . Here 
Z corresponds to the single particle residue and A sp to the 
single particle gap. (b) Size dependence and extrapolation of 
the single particle gap. 



B. Magnetization from pinning fields. 

We have used the pinning field approach to compute 
the staggered moment as a function of U/t. The results at 
ho = 5i are reported in Fig. [3j The extrapolation to the 
thermodynamic limit is carried out using a polynomial 
scaling up to second order in 1/L. Fig. [4] plots the so 
obtained staggered moment for two choices of the pinning 
field as well as the single particle gap. Several comments 
are in order. 

• Within our accuracy, and maybe most importantly, 
with the polynomial fit used in extrapolating the 
data to the thermodynamic limit, it appears that 
the single particle gap opens right when magnetic 
ordering sets in. The only mismatch is at U/t = 3.7 
where we do not detect magnetic ordering but we 




FIG. 3: Magnetic moment at pinning field ho — 5t as a 
function of U/t. Here we have used 9t = 320 and Art = 0.1. 
The solid lines are least square fits to the form a + b/L + c/L . 
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FIG. 4: Staggered moment extrapolated to the thermody- 
namic limit (see Fig. [3J) for two values of the pinning field. 
We have equally plotted the single particle gap in units of U. 
The inset plots the staggered magnetization as obtained from 
a mean-field spin density wave Ansatz. 



do detect a small single particle gap. 

The QMC data in Fig. [4] shows that over a wide pa- 
rameter range, the single particle gap measured in 
units of the Hubbard U, tracks the staggered mag- 
netization. We take this as a strong indication, that 
the magnetization provides the only relevant scale 
in the problem, determining directly the single par- 
ticle gap. We will see below, that this conclusion, 
based here on a simple, polynomial extrapolation 
of the finite size data, is also obtained, if a more 
refined data analysis is performed. 

The data in Fig. [4] exhibits an unusual inflection 
point at approximately U/t = 4.1. Such an inflec- 
tion point is clearly absent at the mean-field level 
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FIG. 5: Data collapse for the magnetization presented in 
Fig. [3] The exponents are taken for the e-expansion of Ref. 
7. (a) The crossing point pins down the value of U c . (b) The 
data collapse, using U c /t — 3.78. 



FIG. 6: Data collapse for the single particle gap. The expo- 
nents are taken for the e-expansion of Ref. ,7. (a) The crossing 
point pins down the value of U c - (b) The data collapse again 
using Uc/t = 3.78. For comparison we have included the data 
for the magnetization. 



(see inset of Fig. Kb. We will discuss the impli- 
cations of this inflection point in the next section. 
Let us finally note, that in previous calculations [1] 
we were unable to resolve staggered moments lesser 
than m ~ 0.03. We thereby missed this inflection 
point in the polynomially extrapolated magnetiza- 
tion curve and concluded the presence of an inter- 
mediate phase. 



C. Finite size scaling 

As mentioned above, one of the particularities of the 
data presented in Fig. [4] is the occurrence of an inflec- 
tion point at U/t = 4.1. It is a natural question to ask 
if this rather peculiar feature may be an artifact of using 
a simple polynomial fitting procedure, which one would 
indeed expect to fail close to criticality. This could result 
in an overestimation of the magnetization in the vicinity 
of the critical point between the semi-metallic and the 
insulating phase of the Hubbard model. As we explain 



next, arguments in favor of this conjecture are provided 
by the large-N treatment of the Gross-Neveu model [5] , 
and the e-expansion around three spatial dimensions in 
the equivalent Gross-Neveu- Yukawa field theory, formu- 
lated in Ref. [7J Given the order parameter exponent, /3, 
as well as the correlation length exponent, v, the stag- 
gered magnetization scales as 



\u-u c f~c 0/v - 



(9) 



Using the standard scaling laws [13], the exponent (3/v 
may conveniently be expressed in terms of the anomalous 
dimension for the order parameter 77, as 



p 



1 



([d+z]-2 + rj). 



(10) 



where d + z is the effective dimensionality of the system. 
If we assume that the Lorentz invariance is emergent at 
the critical point, as it indeed is close to the upper criti- 
cal dimension d up = 3 of the Gross-Neveu- Yukawa theory 
[7], and maybe even more generally |14] , the dynamical 



critical exponent is z = 1. If then the anomalous dimen- 
sion of the order parameter is such that r\ < 3 — d, we find 
that the combination of the exponents [3/v < 1, and our 
polynomial fitting procedure in the previous section could 
very well overestimate the value of the staggered magne- 
tization. In fact both the large-N approach [8] and the 
expansion around the upper critical dimension [7J sug- 
gest that this is indeed the case. Within the first order 
of the expansion in the parameter e = 3 — d, for example, 
■q — 4e/5, so that 
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10 



+ 0(e 2 ). 



(11) 



In two dimensions then, [3/v ~ 0.9. 

To look for the signs of the Gross-Neveu criticality in 
the Hubbard model we have carried out a finite size scal- 
ing analysis based on the usual scaling form 

m = L- p l"F{L 1 / v {U ~U C )). (12) 

Figure km a) plots mJ-P' v versus U for the magnetization 
data at the fixed field ho = 5. (We will omit at this point 
the second scaling variable, hoL y ~ d , since the scaling di- 
mension y — d = (e — rf)/2 <C 1. This second argument 
of the scaling function, present in principle, is therefore 
effectively constant at a fixed ho, and its inclusion does 
not visibly affect the quality of scaling. For further dis- 
cussion of this point, see the Appendix.) As a guide, we 
have used the first-order e-expansion value of (3/v = 0.9. 
Five curves then all cross at a single point U c /t ~ 3.8, 
thereby providing a first non-trivial indication of the crit- 
ical point. This value of U c is slightly larger than that 
obtained with the polynomial fit. The e-expansion value 
of the correlation length exponent reads [7J , 
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(13) 



With this value of v, again at e = 1, we obtain an excel- 
lent data collapse, as shown in Fig. N5[b). 

In accord with the Gross-Neveu- Yukawa theory, the 
numerical data of Fig. [4] support the interpretation that 
the magnetization is the only scale in the problem. To 
further check this interpretation, we have scaled the sin- 
gle single particle gap to the form: 
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*p =L-^ v F{L y ' v {U-U c )) 



(14) 



again using the e-expansion exponents in two dimensions. 
Fig. 6[a) shows that the crossing point of the —ff-L^l v 
curves again occur at U c /t ~ 3.8 and Fig. [6Vb) shows the 
collapse. It is quite remarkable, that within our precision, 
the two scaling functions are equal, F = F. 

Hence, the scaling analysis of our QMC within the 
Gross-Neveu scenario is consistent with a single contin- 
uous quantum phase transition between the semimetal 
and the antiferromagnetic insulator, and suggests that 
the first-order expansion around the upper critical di- 
mension in the Gross-Neveu- Yukawa theory may already 
yield rather accurate values of the critical exponents. 



IV. DISCUSSION AND CONCLUSIONS 

We have introduced an alternative method to com- 
pute the staggered magnetization. By introducing a local 
magnetic field we pin the quantization axis of the ordered 
state. The staggered magnetic moment then corresponds 
to the local magnetization infinitely far away from the 
pinning center. The approach has the major advantage 
that we compute directly the magnetization as opposed 
to its square when measuring correlation functions. We 
can therefore expect improved resolution when the local 
moment is small. One advantage of the approach is an 
internal cross-check which requires the staggered moment 
to be independent on the numerical value of the pinning 
field. We have been able to reach this internal cross-check 
only in the case of large pinning fields. This is consistent 
to the approach proposed by ^ where the pinning field 
is set to infinity on the boundary of the lattice. If the 
pinning field is too small, the approach suffers from large 
and non-monotonous size effects, since the energy scale 
set by the local magnetic field is unable to overcome the 
finite size spin gap. Under these circumstances the ex- 
trapolated value of the magnetization has the tendency 
of underestimating the order. In this article we have 
only considered a very specific form of the pinning field. 
The number of different choices of pinning fields provides 
a playground for optimization of the approach, and for 
minimization of the size effects. 

The application to the Hubbard model on the honey- 
comb lattice sheds new light on the phase diagram of this 
well known problem. The enhanced precision in compar- 
ison to Refs. [T] and [U reveals that the staggered moment 
has the same functional form as the single particle gap 
(as measured in units of U). Remarkably, an excellent 
data collapse onto a single universal curve is found in the 
finite size scaling of both quantities, with the values of 
the critical exponents characteristic of the Gross-Neveu 
criticality between the semimetallic and the magnetic in- 
sulating phases. 
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V. APPENDIX 

Strictly speaking, the finite-size scaling form for the 
magnetization in presence of the external held is 

m = L-P/"G(j,h^) = L-^G(|, !£?), (15) 

where G(X, Y) is the scaling function of two variables, 
and £ is the (diverging) correlation length. Here we have 
used the fact that the relevant Fourier component of the 
local pinning field scales as fy . The dimension y of the 
uniform external field is 1131 



y 



d + z + 2-r] 



(16) 



Away from criticality, £ is bounded and in the thermo- 
dynamic limit, the magnetization scales to its field in- 
dependent value. The data in Fig. 1 confirms this and 
shows that the extrapolated value of the magnetization 
is ho independent provided that, for the considered lat- 
tice sizes, ho is chosen to be large enough. It is also 
interesting to point out that for values of ho in the range 
It < ho < 5t the finite size value of the magnetization is 
next to independent on the value of the pinning field. 
At criticality we can replace £ by L to obtain: 



m = L-^ v G(l,hoL v ~ d ). 



(17) 



With the Lorentz symmetry at the critical point the value 
of the dynamical critical exponent is pinned to z = 1, and 
the scaling dimension of the local field h is 



(18) 



In the e-expansion then, 



y-d=- + 0(e*), 



(19) 



and rather small, presumably even for e = 1 (d = 2). 
For a fixed value of ho, therefore, the second argument 
of the scaling function G is almost constant. To be more 
precise, we can use the asymptotic form (7(1, Y) ex Y 1 / 8 
\ such that in the large- L limit, m ex 



with § = -2- 

6 vy 



Within the 



(y-d) 



0.05 



t u - » . wiuiiiii uiic e expansion — j- 
which results in a very small correction for considered 
lattice sizes. 

Hence on the whole, the scaling function G depends 
rather weakly on the second argument and for practical 
purposes it suffices to neglect it. 
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